float nextIntersection(float4 origin, float4 dir, float4 A, float boxSize);
int packBrickIndex(float x, float y, float z, int borderSize, int brickSize);
float4 getBrickTextureCoords(int bi, int s, int4 ts);

kernel void castRay(global const int* nodePool,
		global const int* dataPool,
		read_only image3d_t brickPool,
		constant const float* vectors,
		constant const float* renderParams,
		global int* image,
		global float* debug) {

	int lid = get_local_id(0);

	// assign the array vectors to cl vector types
	local float4 lfwd, lleft, ldown, origin, nearPlaneCenter;
	local float fov, brickSize, volSize, nearPlaneSize, nearPlaneDist, width, height, pixelSize, lod;
	local int rootNode;
	local int lNodePool[2048], lDataPool[2048];
	local int4 brickPoolDims;
	if (lid == 0) {
		lfwd = vload4(0, vectors);
		lleft = vload4(1, vectors);
		ldown = vload4(2, vectors);
		origin = vload4(3, vectors);
	} else if (lid == 1) {

		fov = renderParams[0];
		brickSize = renderParams[1];
		volSize = renderParams[2];
		nearPlaneSize = renderParams[3];
		nearPlaneDist = renderParams[4];
		width = renderParams[5];
		height = renderParams[6];

		pixelSize = native_divide(nearPlaneSize, width);
		lod = pixelSize * fov * 1.0f;
		nearPlaneCenter = origin + (lfwd * nearPlaneDist);

		rootNode = nodePool[0];
		brickPoolDims = get_image_dim(brickPool);
	}
	barrier(CLK_LOCAL_MEM_FENCE);

	const sampler_t brickSampler = CLK_FILTER_NEAREST | CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE;

	int id = get_global_id(0);

	int2 indices;
	indices.x = id % (int)width;
	indices.y = (int)(id / width);

	// find this pixel on the near plane
	float4 pix = nearPlaneCenter;

	float ri, rj; // i and j relative to the center of the near plane
	ri = indices.x - (width * 0.5f);
	rj = indices.y - (height * 0.5f);

	float di, dj; // pixel's distance from the center
	di = native_divide(ri, width) * nearPlaneSize * 0.5f;
	dj = native_divide(rj, height) * ((nearPlaneSize * native_divide(height, width)) * 0.5f);

	pix -= (di * lleft);
	pix -= (dj * ldown);

	float4 dir = (pix - origin);
	dir = fast_normalize(dir);

	// initialize accumulators

	float4 tclr = 0.0f; // total ray color
	float4 nclr = 0.0f; // n-th sample color

	float totalt = 0.0f; // ray length

	// 4. find the sample position
	float isect;

	float samplet = 1.0f;
	float4 p = pix;

	// setup for octree traversal
	int currentNode; // index of current node
	float4 c; // position offsets in the current node
	float4 corner; // node's lower bound corner
	float nodeSize;  // current node size

	float vSize; // current node voxel size
	float vProjSize; // projected voxel size for this lod
	int currentNodeValue;

	bool searching;

	while ((tclr.s0 < 1.0f) && (samplet > 0.0f)){ // while there are more samples to take
		// 4. set the next sample position
		p += dir * samplet;
		totalt += samplet;
		samplet = 0.0f;

		// descend octree
		currentNode = 0;
		corner = 0.0f;
		c = native_divide(p, volSize); // normalize coordinates
		currentNodeValue = rootNode;

		if (any(c < -0.1f || c > 1.1f)) {
			samplet = nextIntersection(p, dir, corner, volSize + 1);
			continue;
		}

		nodeSize = volSize;

		vSize = native_divide(nodeSize, brickSize);
		vProjSize = vSize * (nearPlaneSize / totalt);

		searching = true;
		// descend until we reach the maximum subdivision level
		// or we reach a node whose size is less than the cone footprint
		while (searching && (vProjSize > lod)) {


			// is this node maximally subdivided?
			if ((currentNodeValue >> 31) || ((currentNodeValue & 0x3FFFFFFF) == 0)) {
				searching = false;
			} else {

				// compute the child node tile offset
				float4 off;
				off = trunc(2.0f * c);

				currentNode = currentNodeValue & 0x3FFFFFFF; // child node tile index
				currentNode += off.x + 2.0f * off.y + 4.0f * off.z;

				c = (2.0f * c) - off;

				nodeSize = native_divide(nodeSize, 2.0f);

				corner += off * (nodeSize);

				vSize = native_divide(nodeSize, brickSize);
				vProjSize = vSize * (nearPlaneSize / totalt);

				currentNodeValue = nodePool[currentNode];	// GLOBAL MEM ACCESS

			}


		}

		// discretized volume rendering integral:
		// C_dst <-- C_dst + (1 - a_dst) * C_src
		// a_dst <-- a_dst + (1 - a_dst) * a_src
		// where opacity a = 1 - T (transparency)

		// find the current node intersection
		isect = nextIntersection(p, dir, corner, nodeSize);


		float step = native_divide(vSize, 1);
		nclr = 0;
		// if this node is not a brick...
		if (!((currentNodeValue >> 30) & 0x1)) {
			// get the node color
			int cl = dataPool[currentNode];		// GLOBAL MEM ACCESS

			float loi = isect;

			/// unpack and correct the color values


			float a = native_divide(cl & 0xFF, 255.0f);			// opacity a
			nclr.s3 = 1.0f - native_powr(1.0f - a, loi);
//			nclr.s3 = a;
			nclr.s0 = nclr.s3 * native_divide((cl >> 8) & 0xFF, 255.0f);	// blue b
			nclr.s1 = nclr.s3 * native_divide((cl >> 16) & 0xFF, 255.0f);	// green g
			nclr.s2 = nclr.s3 * native_divide((cl >> 24) & 0xFF, 255.0f);	// red r

			tclr.s0 += (1.0f - tclr.s3) * nclr.s0;
			tclr.s1 += (1.0f - tclr.s3) * nclr.s1;
			tclr.s2 += (1.0f - tclr.s3) * nclr.s2;

			tclr.s3 += (1.0f - tclr.s3) * nclr.s3;

			// distance to next sample
			samplet += loi + 0.5f;

		} else {

			// find the voxel index
			int poolIndex = dataPool[currentNode];		// GLOBAL MEM ACCESS
			float4 t = getBrickTextureCoords(poolIndex, brickSize, brickPoolDims);

			t += p - corner;

			uint4 vclr;

			while ((tclr.s0 < 1.0f) && (samplet < isect)) {

				vclr = read_imageui(brickPool, brickSampler, t);

				float a = native_divide(vclr.s0 & 0xFF, 255.0f);
				nclr.s3 = a;
				nclr.s0 = a * native_divide(vclr.s1 & 0xFF, 255.0f);
				nclr.s1 = a * native_divide(vclr.s2 & 0xFF, 255.0f);
				nclr.s2 = a * native_divide(vclr.s3 & 0xFF, 255.0f);

				tclr.s0 += (1.0f - tclr.s3) * nclr.s0;
				tclr.s1 += (1.0f - tclr.s3) * nclr.s1;
				tclr.s2 += (1.0f - tclr.s3) * nclr.s2;

				tclr.s3 += (1.0f - tclr.s3) * nclr.s3;

				t += (dir * step);
				samplet += step;
			}

		}

		samplet = (samplet > 0.0f) ? samplet + 0.5f : samplet;

	}

	tclr = clamp(tclr, 0.0f, 1.0f);

//	vstore4(tclr, id, image);

//	write_imagef(image, indices, tclr);

	image[id] = ((int)(tclr.s3 * 255.0f) & 0xFF) |
			(((int)(tclr.s2 * 255.0f) & 0xFF) << 8) |
			(((int)(tclr.s1 * 255.0f) & 0xFF) << 16) |
			(((int)(tclr.s0 * 255.0f) & 0xFF) << 24);

}


float nextIntersection(float4 origin, float4 dir, float4 A, float boxSize) {
	float4 B = A + boxSize;
	float4 invdir = native_recip(dir);
	float isect = -1.0f;

	float tx1 = (A.x - origin.x) * invdir.x;
	float tx2 = (B.x - origin.x) * invdir.x;

	float tmin = min(tx1, tx2);
	float tmax = max(tx1, tx2);

	int cmp = (tmin > tmax);

	float ty1 = (A.y - origin.y) * invdir.y;
	float ty2 = (B.y - origin.y) * invdir.y;

	tmin = max(tmin, min(ty1, ty2));
	tmax = min(tmax, max(ty1, ty2));

	cmp = select(tmin > tmax, cmp, cmp);

	float tz1 = (A.z - origin.z) * invdir.z;
	float tz2 = (B.z - origin.z) * invdir.z;

	tmin = max(tmin, min(tz1, tz2));
	tmax = min(tmax, max(tz1, tz2));

	cmp = select(tmin > tmax, cmp, cmp);
	int cmp2 = (tmin > 0.0f);

	isect = select(tmax, tmin, cmp2);
	isect = select(isect, -1.0f, cmp);

	return isect;
}

float4 getBrickTextureCoords(int bi, int s, int4 ts) {

	int4 b;
	b.x = bi % s;
	b.y = (bi / s) % s;
	b.z = (bi / (s * s)) % s;
	b.w = (bi / (s * s * s));

	int4 t;
	t.x = (b.x + (b.w * s) % ts.x);
	t.y = (b.y + ((b.w / (ts.x / s)) * s) % ts.y);
	t.z = (b.z + ((b.w / ((ts.x / s) * (ts.y / s))) * s) % ts.z);
	t.w = 0.0f;

	return convert_float4(t);
}

int packBrickIndex(float x, float y, float z, int borderSize, int brickSize) {
	clamp(x, 0.0f, brickSize - 1.0f);
	clamp(y, 0.0f, brickSize - 1.0f);
	clamp(z, 0.0f, brickSize - 1.0f);
	int bi = (int)(x + borderSize);
	bi += (int)((y + borderSize) * brickSize);
	bi += (int)((z + borderSize) * brickSize * brickSize);

	return bi;
}
